knitr::opts_chunk$set(warning = FALSE, message = FALSE)
library(spdep) # install.packages("spdep")
Loading required package: sp
Loading required package: spData
To access larger datasets in this package, install the
spDataLarge package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
Loading required package: sf
Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(tidyverse)
Registered S3 method overwritten by 'dplyr':
method from
print.rowwise_df
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
[30m-- [1mAttaching packages[22m --------------------------------------- tidyverse 1.3.0 --[39m
[30m[32mv[30m [34mggplot2[30m 3.3.0 [32mv[30m [34mpurrr [30m 0.3.4
[32mv[30m [34mtibble [30m 3.0.1 [32mv[30m [34mdplyr [30m 0.8.5
[32mv[30m [34mtidyr [30m 1.0.2 [32mv[30m [34mstringr[30m 1.4.0
[32mv[30m [34mreadr [30m 1.3.1 [32mv[30m [34mforcats[30m 0.5.0[39m
[30m-- [1mConflicts[22m ------------------------------------------ tidyverse_conflicts() --
[31mx[30m [34mdplyr[30m::[32mfilter()[30m masks [34mstats[30m::filter()
[31mx[30m [34mdplyr[30m::[32mlag()[30m masks [34mstats[30m::lag()[39m
library(summarytools)
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'pryr':
method from
print.bytes Rcpp
For best results, restart R session and update pander using devtools:: or remotes::install_github('rapporter/pander')
Attaching package: 㤼㸱summarytools㤼㸲
The following object is masked from 㤼㸱package:tibble㤼㸲:
view
library(xtable)
Attaching package: 㤼㸱xtable㤼㸲
The following objects are masked from 㤼㸱package:summarytools㤼㸲:
label, label<-
library(knitr)
library(ExPanDaR)
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'data.table':
method from
print.data.table
library(plotly)
Attaching package: 㤼㸱plotly㤼㸲
The following object is masked from 㤼㸱package:ggplot2㤼㸲:
last_plot
The following object is masked from 㤼㸱package:stats㤼㸲:
filter
The following object is masked from 㤼㸱package:graphics㤼㸲:
layout
library(REAT)
Attaching package: 㤼㸱REAT㤼㸲
The following object is masked from 㤼㸱package:readr㤼㸲:
spec
library(hdrcde)
This is hdrcde 3.3
library(pdfCluster)
pdfCluster 1.0-3
Attaching package: 㤼㸱pdfCluster㤼㸲
The following object is masked from 㤼㸱package:plotly㤼㸲:
groups
The following object is masked from 㤼㸱package:dplyr㤼㸲:
groups
library(readxl)
library(statip)
Attaching package: 㤼㸱statip㤼㸲
The following object is masked from 㤼㸱package:pdfCluster㤼㸲:
plot
The following object is masked from 㤼㸱package:REAT㤼㸲:
cv
The following object is masked from 㤼㸱package:sp㤼㸲:
plot
options(prompt="R> ", digits=3, scipen=999)
library(sf)
cpi <- st_read("Shape file CPI/82 cpi foodstuffs point.shp")
Reading layer `82 cpi foodstuffs point' from data source `D:\Mygithub\tutoring_Felipe\Harry\Shape file CPI\82 cpi foodstuffs point.shp' using driver `ESRI Shapefile'
Simple feature collection with 82 features and 154 fields
geometry type: POINT
dimension: XY
bbox: xmin: 95.3 ymin: -10.2 xmax: 141 ymax: 5.56
geographic CRS: WGS 84
Warning message:
In readChar(file, size, TRUE) : truncating string with embedded nuls
#cpi <- st_read("Shape file CPI/82 cpi procfood point.shp")
#cpi <- st_read("Shape file CPI/82 cpi housing point.shp")
#cpi <- st_read("Shape file CPI/82 cpi health point.shp")
#cpi <- st_read("Shape file CPI/82 cpi clothing point.shp")
#cpi <- st_read("Shape file CPI/82 cpi education point.shp")
#cpi <- st_read("Shape file CPI/82 cpi transport point.shp")
#cpi <- st_read("Shape file CPI/82 cpi general.shp")
plot(cpi)
plotting the first 10 out of 154 attributes; use max.plot = 154 to plot all
class(cpi)
[1] "sf" "data.frame"
head(cpi)
Simple feature collection with 6 features and 154 fields
geometry type: POINT
dimension: XY
bbox: xmin: 113 ymin: -8.69 xmax: 116 ymax: -7.99
geographic CRS: WGS 84
X.x..R Province Regency_Ci Region X Y
1 <NA> Bali Capital City Java -Bali 115 -8.69
2 <NA> West Nusa Tenggara Capital City Nusa Tenggara 116 -8.60
3 <NA> East Java Regency Java -Bali 114 -8.21
4 <NA> East Java Regency Java -Bali 114 -8.16
5 <NA> Bali Regency Java -Bali 115 -8.12
6 <NA> East Java City Java -Bali 113 -7.99
reference IDCity id city products mo2014m1 mo2014m2
1 Kota Denpasar 5171 5171 Denpasar General 109 110
2 Kota Mataram 5271 5271 Mataram General 111 112
3 Banyuwangi 3510 3510 Banyuwangi General 111 112
4 Jember 3509 3509 Jember General 111 111
5 Buleleng 5108 5108 Singaraja General 115 115
6 Kota Malang 3573 3573 Malang General 111 111
mo2014m3 mo2014m4 mo2014m5 mo2014m6 mo2014m7 mo2014m8 mo2014m9
1 110 110 110 110 111 111 112
2 111 111 111 111 112 113 113
3 112 112 112 113 113 113 113
4 111 111 111 111 112 112 112
5 115 115 117 116 117 118 119
6 112 112 112 112 113 114 114
mo2014m10 mo2014m11 mo2014m12 mo2015m1 mo2015m2 mo2015m3 mo2015m4
1 112 114 116 116 116 116 117
2 114 115 117 118 117 118 118
3 113 115 118 118 117 117 117
4 112 114 118 117 117 117 117
5 120 122 125 125 125 126 126
6 114 116 119 119 119 119 120
mo2015m5 mo2015m6 mo2015m7 mo2015m8 mo2015m9 mo2015m10 mo2015m11
1 117 117 119 119 119 118 118
2 118 118 119 119 120 120 120
3 118 118 119 119 119 119 119
4 117 118 119 119 120 119 120
5 127 126 128 128 128 127 127
6 120 121 121 122 122 122 122
mo2015m12 mo2016m1 mo2016m2 mo2016m3 mo2016m4 mo2016m5 mo2016m6
1 120 120 120 120 120 120 121
2 121 123 122 122 122 122 123
3 120 121 121 121 120 121 121
4 120 121 121 121 120 121 121
5 129 131 130 131 131 131 131
6 123 124 124 124 123 123 124
mo2016m7 mo2016m8 mo2016m9 mo2016m10 mo2016m11 mo2016m12 mo2017m1
1 121 122 122 122 122 123 125
2 124 123 123 123 123 124 126
3 122 122 122 122 122 122 123
4 121 121 121 121 121 123 124
5 132 134 134 133 134 135 138
6 125 125 125 125 126 126 128
mo2017m2 mo2017m3 mo2017m4 mo2017m5 mo2017m6 mo2017m7 mo2017m8
1 125 125 125 126 126 126 126
2 127 126 126 126 127 128 127
3 124 123 124 124 125 125 125
4 125 124 125 125 126 126 126
5 139 138 137 137 136 137 137
6 128 128 129 130 130 131 130
mo2017m9 mo2017m10 mo2017m11 mo2017m12 mo2018m1 mo2018m2 mo2018m3
1 126 126 126 127 128 129 129
2 127 128 128 129 129 130 130
3 125 125 126 126 127 127 128
4 126 126 126 127 128 128 128
5 136 136 138 140 141 141 142
6 130 130 130 131 132 132 132
mo2018m4 mo2018m5 mo2018m6 mo2018m7 mo2018m8 mo2018m9 mo2018m10
1 129 129 130 131 131 130 130
2 130 130 131 132 132 131 132
3 128 128 128 129 128 128 128
4 128 129 130 129 129 129 130
5 141 141 141 142 142 141 141
6 133 133 133 134 134 133 134
mo2018m11 mo2018m12 mo2019m1 mo2019m2 mo2019m3 mo2019m4 mo2019m5
1 130 132 132 132 132 132 133
2 132 133 133 133 133 133 134
3 128 129 129 129 130 130 130
4 130 131 131 131 131 131 132
5 141 142 143 143 143 144 144
6 134 135 136 135 136 136 137
mo2019m6 mo2019m7 mo2019m8 mo2019m9 mo2019m10 mo2019m11 mo2019m12
1 133 134 134 133 134 134 135
2 135 135 135 134 135 135 135
3 131 131 131 131 131 132 132
4 132 132 132 132 132 133 133
5 144 146 146 145 145 145 146
6 136 137 137 137 137 137 138
ln2014m1 ln2014m2 ln2014m3 ln2014m4 ln2014m5 ln2014m6 ln2014m7
1 4.69 4.70 4.70 4.70 4.70 4.70 4.71
2 4.71 4.71 4.71 4.71 4.71 4.71 4.72
3 4.71 4.72 4.72 4.72 4.72 4.72 4.73
4 4.71 4.71 4.71 4.71 4.71 4.71 4.72
5 4.74 4.75 4.75 4.75 4.76 4.76 4.76
6 4.71 4.71 4.72 4.72 4.72 4.72 4.73
ln2014m8 ln2014m9 ln2014m10 ln2014m11 ln2014m12 ln2015m1 ln2015m2
1 4.71 4.72 4.72 4.74 4.76 4.76 4.76
2 4.73 4.73 4.73 4.74 4.77 4.77 4.77
3 4.72 4.73 4.73 4.74 4.77 4.77 4.76
4 4.72 4.72 4.72 4.74 4.77 4.76 4.76
5 4.77 4.78 4.78 4.80 4.83 4.83 4.83
6 4.73 4.73 4.74 4.75 4.78 4.78 4.78
ln2015m3 ln2015m4 ln2015m5 ln2015m6 ln2015m7 ln2015m8 ln2015m9
1 4.76 4.76 4.76 4.77 4.77 4.78 4.78
2 4.77 4.77 4.77 4.77 4.78 4.78 4.79
3 4.76 4.76 4.77 4.77 4.78 4.78 4.78
4 4.76 4.76 4.77 4.77 4.78 4.78 4.78
5 4.83 4.84 4.84 4.84 4.85 4.85 4.85
6 4.78 4.78 4.79 4.79 4.80 4.80 4.80
ln2015m10 ln2015m11 ln2015m12 ln2016m1 ln2016m2 ln2016m3 ln2016m4
1 4.77 4.77 4.78 4.79 4.79 4.79 4.79
2 4.79 4.79 4.80 4.81 4.81 4.81 4.80
3 4.78 4.78 4.79 4.80 4.80 4.80 4.79
4 4.78 4.79 4.79 4.79 4.80 4.80 4.79
5 4.84 4.85 4.86 4.87 4.87 4.88 4.88
6 4.80 4.80 4.81 4.82 4.82 4.82 4.81
ln2016m5 ln2016m6 ln2016m7 ln2016m8 ln2016m9 ln2016m10 ln2016m11
1 4.79 4.79 4.80 4.80 4.81 4.80 4.81
2 4.80 4.81 4.82 4.82 4.81 4.81 4.82
3 4.79 4.80 4.80 4.80 4.80 4.80 4.80
4 4.79 4.80 4.80 4.80 4.80 4.80 4.80
5 4.88 4.88 4.89 4.89 4.90 4.89 4.90
6 4.82 4.82 4.83 4.83 4.83 4.83 4.83
ln2016m12 ln2017m1 ln2017m2 ln2017m3 ln2017m4 ln2017m5 ln2017m6
1 4.81 4.83 4.83 4.83 4.83 4.83 4.83
2 4.82 4.84 4.84 4.84 4.83 4.84 4.84
3 4.81 4.81 4.82 4.82 4.82 4.82 4.83
4 4.81 4.82 4.83 4.82 4.83 4.83 4.83
5 4.91 4.92 4.93 4.93 4.92 4.92 4.92
6 4.84 4.85 4.86 4.85 4.86 4.87 4.87
ln2017m7 ln2017m8 ln2017m9 ln2017m10 ln2017m11 ln2017m12 ln2018m1
1 4.83 4.84 4.83 4.83 4.83 4.85 4.85
2 4.85 4.85 4.85 4.85 4.85 4.86 4.86
3 4.83 4.83 4.83 4.83 4.83 4.84 4.85
4 4.84 4.83 4.83 4.83 4.84 4.84 4.85
5 4.92 4.92 4.91 4.91 4.93 4.94 4.95
6 4.87 4.87 4.87 4.87 4.87 4.88 4.88
ln2018m2 ln2018m3 ln2018m4 ln2018m5 ln2018m6 ln2018m7 ln2018m8
1 4.86 4.86 4.86 4.86 4.87 4.87 4.87
2 4.87 4.87 4.87 4.87 4.87 4.88 4.88
3 4.85 4.85 4.85 4.85 4.86 4.86 4.86
4 4.85 4.85 4.85 4.86 4.86 4.86 4.86
5 4.95 4.95 4.95 4.95 4.95 4.95 4.96
6 4.88 4.89 4.89 4.89 4.89 4.89 4.90
ln2018m9 ln2018m10 ln2018m11 ln2018m12 ln2019m1 ln2019m2 ln2019m3
1 4.87 4.87 4.87 4.88 4.89 4.88 4.88
2 4.88 4.88 4.88 4.89 4.89 4.89 4.89
3 4.85 4.85 4.85 4.86 4.86 4.86 4.86
4 4.86 4.86 4.87 4.87 4.87 4.87 4.87
5 4.95 4.95 4.95 4.96 4.96 4.96 4.96
6 4.89 4.89 4.90 4.91 4.91 4.91 4.91
ln2019m4 ln2019m5 ln2019m6 ln2019m7 ln2019m8 ln2019m9 ln2019m10
1 4.89 4.89 4.89 4.89 4.90 4.89 4.89
2 4.89 4.90 4.91 4.91 4.90 4.90 4.90
3 4.87 4.87 4.87 4.88 4.88 4.88 4.88
4 4.88 4.88 4.88 4.88 4.89 4.88 4.88
5 4.97 4.97 4.97 4.98 4.99 4.98 4.98
6 4.91 4.92 4.92 4.92 4.92 4.92 4.92
ln2019m11 geometry
1 4.89 POINT (115 -8.69)
2 4.90 POINT (116 -8.6)
3 4.88 POINT (114 -8.21)
4 4.89 POINT (114 -8.16)
5 4.98 POINT (115 -8.12)
6 4.92 POINT (113 -7.99)
tail(cpi)
Simple feature collection with 6 features and 154 fields
geometry type: POINT
dimension: XY
bbox: xmin: 120 ymin: -0.919 xmax: 134 ymax: 1.56
geographic CRS: WGS 84
X.x..R Province Regency_Ci Region X Y
77 <NA> Papua Barat Capital City Papua 131 -0.919
78 <NA> Papua Barat Regency Papua 134 -0.860
79 <NA> Central Sulawesi Capital City Sulawesi 120 -0.808
80 <NA> Gorontalo Capital City Sulawesi 123 0.547
81 <NA> North Maluku Capital City Maluku 127 0.811
82 <NA> North Sulawesi Capital City Sulawesi 125 1.563
reference IDCity id city products mo2014m1
77 Kota Sorong 9171 9171 Sorong General 108
78 Manokwari 9105 9105 Manokwari General 106
79 Kota Palu 7271 7271 Palu General 112
80 Kota Gorontalo, Kota 7571 7571 Gorontalo General 109
81 Kota Ternate, Kota 8271 8271 Ternate General 112
82 Kota Manado 7171 7171 Manado General 109
mo2014m2 mo2014m3 mo2014m4 mo2014m5 mo2014m6 mo2014m7 mo2014m8
77 109 109 110 110 110 112 114
78 107 106 106 107 107 108 110
79 111 111 112 113 114 115 116
80 108 108 109 109 109 110 110
81 112 112 113 113 114 117 116
82 109 109 110 110 110 111 111
mo2014m9 mo2014m10 mo2014m11 mo2014m12 mo2015m1 mo2015m2 mo2015m3
77 115 114 114 116 116 117 117
78 110 111 111 113 112 112 113
79 115 117 117 120 120 118 117
80 110 110 111 115 114 113 114
81 117 118 119 122 122 121 121
82 111 112 114 119 118 118 118
mo2015m4 mo2015m5 mo2015m6 mo2015m7 mo2015m8 mo2015m9 mo2015m10
77 117 117 120 122 123 123 123
78 113 113 114 115 113 114 113
79 118 120 120 122 121 121 122
80 114 115 116 117 118 118 118
81 122 123 124 125 127 125 126
82 118 119 120 121 121 121 123
mo2015m11 mo2015m12 mo2016m1 mo2016m2 mo2016m3 mo2016m4 mo2016m5
77 122 123 125 125 125 124 123
78 113 116 116 116 116 116 117
79 123 125 125 124 124 124 125
80 118 120 120 120 120 120 120
81 126 128 128 127 128 128 128
82 123 125 125 124 124 123 123
mo2016m6 mo2016m7 mo2016m8 mo2016m9 mo2016m10 mo2016m11 mo2016m12
77 124 126 127 127 126 126 127
78 119 120 122 121 120 121 122
79 126 126 126 126 125 126 127
80 122 122 121 121 120 121 122
81 128 130 130 130 130 130 130
82 124 125 125 124 124 128 126
mo2017m1 mo2017m2 mo2017m3 mo2017m4 mo2017m5 mo2017m6 mo2017m7
77 128 128 129 128 128 129 130
78 122 122 122 121 122 124 125
79 129 129 129 130 131 132 132
80 123 124 124 124 124 126 127
81 131 131 131 131 131 133 135
82 127 128 129 129 127 129 130
mo2017m8 mo2017m9 mo2017m10 mo2017m11 mo2017m12 mo2018m1 mo2018m2
77 129 129 129 128 129 129 130
78 123 125 124 124 125 126 124
79 132 132 130 130 133 134 133
80 126 126 126 126 127 128 127
81 133 132 133 131 133 134 134
82 130 128 128 128 129 129 130
mo2018m3 mo2018m4 mo2018m5 mo2018m6 mo2018m7 mo2018m8 mo2018m9
77 130 131 132 134 136 136 135
78 125 125 126 127 128 128 128
79 133 134 134 137 137 137 135
80 127 127 128 129 129 129 129
81 135 136 136 139 137 137 137
82 130 132 132 133 132 131 130
mo2018m10 mo2018m11 mo2018m12 mo2019m1 mo2019m2 mo2019m3 mo2019m4
77 135 135 135 135 134 133 134
78 129 130 132 133 133 133 133
79 138 140 141 141 141 140 141
80 129 129 130 130 129 129 130
81 137 137 138 139 139 139 139
82 130 133 134 135 134 133 132
mo2019m5 mo2019m6 mo2019m7 mo2019m8 mo2019m9 mo2019m10 mo2019m11
77 135 136 136 137 137 137 135
78 136 136 135 136 136 137 137
79 143 144 143 144 143 143 143
80 132 132 132 133 133 133 133
81 140 141 141 141 140 140 141
82 135 140 138 136 135 136 141
mo2019m12 ln2014m1 ln2014m2 ln2014m3 ln2014m4 ln2014m5 ln2014m6
77 136 4.69 4.69 4.69 4.70 4.70 4.70
78 138 4.67 4.67 4.67 4.67 4.67 4.68
79 144 4.71 4.71 4.71 4.72 4.72 4.73
80 134 4.69 4.68 4.68 4.69 4.69 4.69
81 141 4.72 4.71 4.72 4.73 4.73 4.74
82 138 4.69 4.69 4.69 4.70 4.70 4.70
ln2014m7 ln2014m8 ln2014m9 ln2014m10 ln2014m11 ln2014m12 ln2015m1
77 4.72 4.74 4.75 4.74 4.74 4.75 4.76
78 4.69 4.70 4.70 4.71 4.71 4.72 4.72
79 4.75 4.75 4.75 4.76 4.76 4.79 4.79
80 4.70 4.70 4.70 4.70 4.71 4.75 4.73
81 4.76 4.75 4.76 4.77 4.78 4.81 4.80
82 4.71 4.71 4.71 4.72 4.74 4.78 4.77
ln2015m2 ln2015m3 ln2015m4 ln2015m5 ln2015m6 ln2015m7 ln2015m8
77 4.76 4.76 4.76 4.77 4.78 4.80 4.81
78 4.72 4.73 4.72 4.72 4.74 4.75 4.73
79 4.77 4.77 4.77 4.79 4.79 4.80 4.80
80 4.73 4.74 4.74 4.75 4.75 4.76 4.77
81 4.79 4.80 4.80 4.81 4.82 4.83 4.84
82 4.77 4.77 4.77 4.78 4.79 4.80 4.79
ln2015m9 ln2015m10 ln2015m11 ln2015m12 ln2016m1 ln2016m2 ln2016m3
77 4.81 4.81 4.81 4.81 4.82 4.83 4.82
78 4.73 4.73 4.73 4.75 4.75 4.75 4.75
79 4.80 4.81 4.81 4.83 4.83 4.82 4.82
80 4.77 4.77 4.77 4.79 4.78 4.79 4.79
81 4.83 4.84 4.84 4.85 4.86 4.85 4.85
82 4.80 4.81 4.81 4.83 4.83 4.82 4.82
ln2016m4 ln2016m5 ln2016m6 ln2016m7 ln2016m8 ln2016m9 ln2016m10
77 4.82 4.81 4.82 4.83 4.85 4.85 4.84
78 4.75 4.76 4.78 4.79 4.80 4.79 4.79
79 4.82 4.83 4.83 4.84 4.83 4.84 4.83
80 4.79 4.79 4.80 4.80 4.80 4.80 4.79
81 4.85 4.85 4.86 4.87 4.86 4.87 4.86
82 4.81 4.81 4.82 4.83 4.83 4.82 4.82
ln2016m11 ln2016m12 ln2017m1 ln2017m2 ln2017m3 ln2017m4 ln2017m5
77 4.84 4.84 4.85 4.85 4.86 4.85 4.85
78 4.80 4.81 4.81 4.80 4.80 4.80 4.81
79 4.83 4.84 4.86 4.86 4.86 4.87 4.88
80 4.80 4.80 4.81 4.82 4.82 4.82 4.82
81 4.87 4.87 4.88 4.88 4.87 4.88 4.88
82 4.85 4.83 4.84 4.86 4.86 4.86 4.85
ln2017m6 ln2017m7 ln2017m8 ln2017m9 ln2017m10 ln2017m11 ln2017m12
77 4.86 4.86 4.86 4.86 4.86 4.85 4.86
78 4.82 4.83 4.81 4.83 4.82 4.82 4.82
79 4.88 4.88 4.88 4.88 4.87 4.87 4.89
80 4.84 4.85 4.84 4.84 4.84 4.84 4.84
81 4.89 4.90 4.89 4.88 4.89 4.88 4.89
82 4.86 4.87 4.86 4.85 4.85 4.85 4.86
ln2018m1 ln2018m2 ln2018m3 ln2018m4 ln2018m5 ln2018m6 ln2018m7
77 4.86 4.87 4.87 4.88 4.88 4.90 4.91
78 4.83 4.82 4.82 4.83 4.84 4.85 4.85
79 4.89 4.89 4.89 4.90 4.90 4.92 4.92
80 4.85 4.84 4.85 4.85 4.85 4.86 4.86
81 4.90 4.90 4.91 4.91 4.92 4.93 4.92
82 4.86 4.87 4.87 4.88 4.89 4.89 4.89
ln2018m8 ln2018m9 ln2018m10 ln2018m11 ln2018m12 ln2019m1 ln2019m2
77 4.92 4.90 4.90 4.91 4.90 4.91 4.90
78 4.85 4.85 4.86 4.87 4.88 4.89 4.89
79 4.92 4.91 4.93 4.94 4.95 4.95 4.95
80 4.86 4.86 4.86 4.86 4.87 4.87 4.86
81 4.92 4.92 4.92 4.92 4.93 4.94 4.93
82 4.88 4.87 4.87 4.89 4.90 4.91 4.90
ln2019m3 ln2019m4 ln2019m5 ln2019m6 ln2019m7 ln2019m8 ln2019m9
77 4.89 4.89 4.91 4.91 4.92 4.92 4.92
78 4.89 4.89 4.91 4.91 4.91 4.92 4.91
79 4.94 4.95 4.96 4.97 4.96 4.97 4.96
80 4.86 4.87 4.88 4.89 4.89 4.89 4.89
81 4.93 4.94 4.94 4.95 4.95 4.95 4.94
82 4.89 4.88 4.91 4.94 4.93 4.91 4.90
ln2019m10 ln2019m11 geometry
77 4.92 4.91 POINT (131 -0.919)
78 4.92 4.92 POINT (134 -0.86)
79 4.96 4.96 POINT (120 -0.808)
80 4.89 4.89 POINT (123 0.547)
81 4.94 4.95 POINT (127 0.811)
82 4.92 4.95 POINT (125 1.56)
#minimum distance is 338 for the dnearneigh fuction
#convert cpi spatial object to data.frame
cpi.df <- as.data.frame(cpi)
nb <- dnearneigh(cpi, d1=0, d2=672+ 10*jj, longlat=T, row.names = IDs) #create weight matrix and name W.matrix
# (a) coordinates and IDs
#coords <- cpi[,c("X", "Y")]
#coords <- coordinates(coords)
IDs <- cpi$IDCity
# (b) identify neighbors given the distance threshold
nb672km <- dnearneigh(cpi, d1=0, d2=672, row.names = IDs)
summary(nb672km)
Neighbour list object:
Number of regions: 82
Number of nonzero links: 1588
Percentage nonzero weights: 23.6
Average number of links: 19.4
Link number distribution:
1 2 3 4 5 6 7 9 10 12 13 14 15 16 17 18 19 20 23 24 25
2 1 1 5 1 1 2 2 2 1 3 2 3 7 1 5 4 2 3 2 3
26 27 28 29 30 31 32 33 34
3 5 3 8 2 3 2 1 2
2 least connected regions:
9401 9471 with 1 link
2 most connected regions:
3374 3319 with 34 links
# (c) calculate the spatial weights matrix
W.matrixdist <- nb2listw(nb672km)
summary(W.matrixdist)
Characteristics of weights list object:
Neighbour list object:
Number of regions: 82
Number of nonzero links: 1588
Percentage nonzero weights: 23.6
Average number of links: 19.4
Link number distribution:
1 2 3 4 5 6 7 9 10 12 13 14 15 16 17 18 19 20 23 24 25
2 1 1 5 1 1 2 2 2 1 3 2 3 7 1 5 4 2 3 2 3
26 27 28 29 30 31 32 33 34
3 5 3 8 2 3 2 1 2
2 least connected regions:
9401 9471 with 1 link
2 most connected regions:
3374 3319 with 34 links
Weights style: W
Weights constants summary:
#nb2listw(nb338km) #, zero.policy = TRUE)
#W.matrix <- nb2listw(nb338km,zero.policy = TRUE)
#summary(W.matrix)
#weight matrix k nearest
# (a) coordinates and IDs
#coords <- cpi[,c("X", "Y")]
#coords <- coordinates(coords)
IDs <- cpi$IDCity
# (b) identify neighbors given the distance threshold
nearneigh4 <- knearneigh(cpi, k=4)
summary(nearneigh4)
Length Class Mode
nn 328 -none- numeric
np 1 -none- numeric
k 1 -none- numeric
dimension 1 -none- numeric
x 164 -none- numeric
# (c) calculate the spatial weights matrix
W.matrixnearest <- nb2listw(knn2nb(nearneigh4))
summary(W.matrixnearest)
Characteristics of weights list object:
Neighbour list object:
Number of regions: 82
Number of nonzero links: 328
Percentage nonzero weights: 4.88
Average number of links: 4
Non-symmetric neighbours list
Link number distribution:
4
82
82 least connected regions:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 with 4 links
82 most connected regions:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 with 4 links
Weights style: W
Weights constants summary:
#checking the average of Moran’s I
cpi <- cpi %>%
mutate(avg = (mo2016m10+mo2016m11+mo2016m12+mo2017m1+mo2017m2+mo2017m3)/6)
cpi.df
cpi[,11]
Simple feature collection with 82 features and 1 field
geometry type: POINT
dimension: XY
bbox: xmin: 95.3 ymin: -10.2 xmax: 141 ymax: 5.56
CRS: 4326
First 10 features:
mo2014m1 geometry
1 116 POINT (115 -8.69)
2 120 POINT (116 -8.6)
3 127 POINT (114 -8.21)
4 117 POINT (114 -8.16)
5 118 POINT (115 -8.12)
6 118 POINT (113 -7.99)
7 116 POINT (112 -7.83)
8 120 POINT (110 -7.79)
9 117 POINT (113 -7.75)
10 120 POINT (109 -7.7)
# (b) identify neighbors given the distance threshold
nb672km <- dnearneigh(cpi, d1=0, d2=672, row.names = IDs)
summary(nb672km)
Neighbour list object:
Number of regions: 82
Number of nonzero links: 1588
Percentage nonzero weights: 23.6
Average number of links: 19.4
Link number distribution:
1 2 3 4 5 6 7 9 10 12 13 14 15 16 17 18 19 20 23 24 25
2 1 1 5 1 1 2 2 2 1 3 2 3 7 1 5 4 2 3 2 3
26 27 28 29 30 31 32 33 34
3 5 3 8 2 3 2 1 2
2 least connected regions:
9401 9471 with 1 link
2 most connected regions:
3374 3319 with 34 links
# (c) calculate the spatial weights matrix
W.matrix <- nb2listw(nb672km)
summary(W.matrix)
Characteristics of weights list object:
Neighbour list object:
Number of regions: 82
Number of nonzero links: 1588
Percentage nonzero weights: 23.6
Average number of links: 19.4
Link number distribution:
1 2 3 4 5 6 7 9 10 12 13 14 15 16 17 18 19 20 23 24 25
2 1 1 5 1 1 2 2 2 1 3 2 3 7 1 5 4 2 3 2 3
26 27 28 29 30 31 32 33 34
3 5 3 8 2 3 2 1 2
2 least connected regions:
9401 9471 with 1 link
2 most connected regions:
3374 3319 with 34 links
Weights style: W
Weights constants summary:
?localG
localG(cpi.df$mo2014m3, W.matrix)
[1] 1.64878 1.38398 1.15990 1.85607 1.38404 2.23471
[7] 2.74276 4.20782 1.88619 4.18377 2.90712 4.10851
[13] 4.11783 2.96376 2.32004 1.54396 3.87443 3.04638
[19] 2.84404 3.55592 3.73188 3.44055 2.83751 2.74521
[25] 2.74096 2.59241 2.64536 2.66724 2.59283 3.15398
[31] 2.79039 2.11122 0.33211 2.41516 2.63926 3.27713
[37] 1.20955 -1.13287 1.12867 -2.04853 -1.28549 -0.92493
[43] -2.56634 -3.12814 -1.45965 -3.13347 -2.42568 -2.84249
[49] -2.48415 -3.21206 -2.68763 0.34423 0.81283 2.67875
[55] -1.31488 2.28088 1.37271 -1.58224 -0.27780 -2.16889
[61] 0.63526 -0.22797 -0.57143 -0.00653 0.28373 -0.05982
[67] 0.06157 -0.27498 -0.73210 -0.32900 -0.30385 -0.64036
[73] -0.35039 -1.28095 -1.35349 -1.37972 -1.34818 -0.83486
[79] -2.45196 -2.64425 -1.16874 -3.28812
attr(,"gstari")
[1] FALSE
attr(,"call")
localG(x = cpi.df$mo2014m3, listw = W.matrix)
attr(,"class")
[1] "localG"
getis<- localG(cpi.df$mo2014m3, nb2listw(nb672km, style="B"), return_internals = T, GeoDa = T)
getis
[1] 1.64878 1.38398 1.15990 1.85607 1.38404 2.23471
[7] 2.74276 4.20782 1.88619 4.18377 2.90712 4.10851
[13] 4.11783 2.96376 2.32004 1.54396 3.87443 3.04638
[19] 2.84404 3.55592 3.73188 3.44055 2.83751 2.74521
[25] 2.74096 2.59241 2.64536 2.66724 2.59283 3.15398
[31] 2.79039 2.11122 0.33211 2.41516 2.63926 3.27713
[37] 1.20955 -1.13287 1.12867 -2.04853 -1.28549 -0.92493
[43] -2.56634 -3.12814 -1.45965 -3.13347 -2.42568 -2.84249
[49] -2.48415 -3.21206 -2.68763 0.34423 0.81283 2.67875
[55] -1.31488 2.28088 1.37271 -1.58224 -0.27780 -2.16889
[61] 0.63526 -0.22797 -0.57143 -0.00653 0.28373 -0.05982
[67] 0.06157 -0.27498 -0.73210 -0.32900 -0.30385 -0.64036
[73] -0.35039 -1.28095 -1.35349 -1.37972 -1.34818 -0.83486
[79] -2.45196 -2.64425 -1.16874 -3.28812
attr(,"internals")
G EG VG
[1,] 0.2261 0.2222 0.000005453
[2,] 0.2503 0.2469 0.000005852
[3,] 0.2868 0.2840 0.000005949
[4,] 0.2887 0.2840 0.000006415
[5,] 0.2875 0.2840 0.000006406
[6,] 0.3392 0.3333 0.000007000
[7,] 0.3654 0.3580 0.000007245
[8,] 0.4066 0.3951 0.000007499
[9,] 0.3011 0.2963 0.000006579
[10,] 0.3693 0.3580 0.000007242
[11,] 0.3907 0.3827 0.000007454
[12,] 0.4187 0.4074 0.000007611
[13,] 0.3815 0.3704 0.000007336
[14,] 0.3660 0.3580 0.000007247
[15,] 0.3767 0.3704 0.000007347
[16,] 0.3126 0.3086 0.000006730
[17,] 0.4305 0.4198 0.000007681
[18,] 0.3538 0.3457 0.000007129
[19,] 0.3533 0.3457 0.000007128
[20,] 0.3924 0.3827 0.000007441
[21,] 0.4297 0.4198 0.000007071
[22,] 0.3921 0.3827 0.000007435
[23,] 0.3657 0.3580 0.000007252
[24,] 0.3654 0.3580 0.000007191
[25,] 0.3654 0.3580 0.000007185
[26,] 0.3402 0.3333 0.000006951
[27,] 0.3527 0.3457 0.000007047
[28,] 0.3279 0.3210 0.000006811
[29,] 0.3402 0.3333 0.000007012
[30,] 0.3417 0.3333 0.000007006
[31,] 0.3400 0.3333 0.000005643
[32,] 0.3017 0.2963 0.000006553
[33,] 0.3219 0.3210 0.000006848
[34,] 0.3149 0.3086 0.000006674
[35,] 0.3651 0.3580 0.000007125
[36,] 0.4040 0.3951 0.000007473
[37,] 0.3118 0.3086 0.000006732
[38,] 0.0480 0.0494 0.000001459
[39,] 0.0130 0.0123 0.000000384
[40,] 0.1076 0.1111 0.000002996
[41,] 0.2192 0.2222 0.000005450
[42,] 0.0361 0.0370 0.000001124
[43,] 0.2162 0.2222 0.000005454
[44,] 0.1540 0.1605 0.000004245
[45,] 0.2311 0.2346 0.000005600
[46,] 0.1905 0.1975 0.000004994
[47,] 0.1557 0.1605 0.000003948
[48,] 0.1912 0.1975 0.000004937
[49,] 0.0584 0.0617 0.000001819
[50,] 0.1782 0.1852 0.000004757
[51,] 0.1916 0.1975 0.000004823
[52,] 0.0126 0.0123 0.000000379
[53,] 0.2365 0.2346 0.000005588
[54,] 0.3652 0.3580 0.000007077
[55,] 0.1946 0.1975 0.000004996
[56,] 0.3270 0.3210 0.000006850
[57,] 0.2378 0.2346 0.000005625
[58,] 0.1572 0.1605 0.000004240
[59,] 0.1969 0.1975 0.000004981
[60,] 0.1682 0.1728 0.000004511
[61,] 0.2484 0.2469 0.000005539
[62,] 0.1970 0.1975 0.000004994
[63,] 0.1101 0.1111 0.000003116
[64,] 0.1975 0.1975 0.000004976
[65,] 0.1240 0.1235 0.000003413
[66,] 0.2221 0.2222 0.000005329
[67,] 0.2347 0.2346 0.000005662
[68,] 0.2092 0.2099 0.000005206
[69,] 0.2205 0.2222 0.000005454
[70,] 0.1845 0.1852 0.000004758
[71,] 0.1722 0.1728 0.000004508
[72,] 0.0486 0.0494 0.000001476
[73,] 0.1474 0.1481 0.000003971
[74,] 0.0844 0.0864 0.000002439
[75,] 0.0843 0.0864 0.000002455
[76,] 0.0721 0.0741 0.000002123
[77,] 0.0478 0.0494 0.000001465
[78,] 0.0240 0.0247 0.000000741
[79,] 0.1798 0.1852 0.000004735
[80,] 0.1187 0.1235 0.000003268
[81,] 0.0480 0.0494 0.000001427
[82,] 0.0454 0.0494 0.000001481
attr(,"gstari")
[1] FALSE
attr(,"call")
localG(x = cpi.df$mo2014m3, listw = nb2listw(nb672km, style = "B"),
return_internals = T, GeoDa = T)
attr(,"class")
[1] "localG"
a<- localG(dist$RGDPPC2017, nb2listw(nb338km, style="B"), return_internals = T, GeoDa = T)
#a
a[1]
#attributes(a)
b<-attributes(a)$internals
#b
mo2014m1.g <- localG(cpi.df$mo2014m1, nb2listw(nb672km, style="B"), return_internals = T, GeoDa = T)
Getis.m2014m1 <- as.data.frame(attr(mo2014m1.g, which = "internals")) # retrieve "internals"
cpi.df$g.2014m1 <- cpi.df$mo2014m1 * (Getis.m2014m1$EG/Getis.m2014m1$G ) # "multiplicative filter"
filtered <- list()
a=c()
a<- cpi.df[,7]
a
[1] 5171 5271 3510 3509 5108 3573 3571 3471 3574 3301 3577
[12] 3372 3302 3278 3578 3529 3374 3272 3273 3376 3319 3274
[23] 3271 3276 3275 3671 3173 3673 3672 1871 1872 1771 6371
[34] 1674 1671 1902 6202 5371 9401 5310 5272 8172 7302 7472
[45] 7371 7311 7471 7372 8171 7373 7604 9471 6271 1971 6309
[56] 1571 1509 6471 1371 6472 1403 1375 6171 1471 6172 2172
[67] 2171 1277 1473 1271 1273 6571 1275 1107 1174 1171 9171
[78] 9105 7271 7571 8271 7171
for (i in 12:82) {
g <- localG(cpi.df[,i], nb2listw(nb672km, style="B"), return_internals = T, GeoDa = T)
Getis <- as.data.frame(attr(g, which = "internals"))
filtered[[i]] <- cpi.df[,i] * (Getis$EG/Getis$G ) # "multiplicative filter"
a<- as.data.frame(cbind(a,filtered[[i]]))
colnames(a)[ncol(a)] <- paste("f",colnames(cpi.df)[i],sep = "_")
}
a
combined_cpi_g <- left_join(cpi.df,a, by = c("IDCity" ="a"))
combined_cpi_g
#saving the filtered data in a new file
dist
dist
out<-dist %>%
select(1:4, 45,46,47, 50, 54, 58, 60,61,62, 68:85)
out
write.csv(out, "01-raw-data/output_spatial_fitlering/2017filt.csv")
s2.df <- data.frame( V1=0, cpigetis= 0)
s2.df
NA
#create a data frame of the same type as the one created by the PPA program # please be aware that the code is much longer and complicated than needed but it was better to use what i have done in the past to study the output of the PPA program
W.matrixnearest <- nb2listw(knn2nb(nearneigh4))
cpigetis <- localG(cpi.df$mo2014m3, nb2listw(nb672km, style="B"), return_internals = T, GeoDa = T)
cpigetis
[1] 1.64878 1.38398 1.15990 1.85607 1.38404 2.23471
[7] 2.74276 4.20782 1.88619 4.18377 2.90712 4.10851
[13] 4.11783 2.96376 2.32004 1.54396 3.87443 3.04638
[19] 2.84404 3.55592 3.73188 3.44055 2.83751 2.74521
[25] 2.74096 2.59241 2.64536 2.66724 2.59283 3.15398
[31] 2.79039 2.11122 0.33211 2.41516 2.63926 3.27713
[37] 1.20955 -1.13287 1.12867 -2.04853 -1.28549 -0.92493
[43] -2.56634 -3.12814 -1.45965 -3.13347 -2.42568 -2.84249
[49] -2.48415 -3.21206 -2.68763 0.34423 0.81283 2.67875
[55] -1.31488 2.28088 1.37271 -1.58224 -0.27780 -2.16889
[61] 0.63526 -0.22797 -0.57143 -0.00653 0.28373 -0.05982
[67] 0.06157 -0.27498 -0.73210 -0.32900 -0.30385 -0.64036
[73] -0.35039 -1.28095 -1.35349 -1.37972 -1.34818 -0.83486
[79] -2.45196 -2.64425 -1.16874 -3.28812
attr(,"internals")
G EG VG
[1,] 0.2261 0.2222 0.000005453
[2,] 0.2503 0.2469 0.000005852
[3,] 0.2868 0.2840 0.000005949
[4,] 0.2887 0.2840 0.000006415
[5,] 0.2875 0.2840 0.000006406
[6,] 0.3392 0.3333 0.000007000
[7,] 0.3654 0.3580 0.000007245
[8,] 0.4066 0.3951 0.000007499
[9,] 0.3011 0.2963 0.000006579
[10,] 0.3693 0.3580 0.000007242
[11,] 0.3907 0.3827 0.000007454
[12,] 0.4187 0.4074 0.000007611
[13,] 0.3815 0.3704 0.000007336
[14,] 0.3660 0.3580 0.000007247
[15,] 0.3767 0.3704 0.000007347
[16,] 0.3126 0.3086 0.000006730
[17,] 0.4305 0.4198 0.000007681
[18,] 0.3538 0.3457 0.000007129
[19,] 0.3533 0.3457 0.000007128
[20,] 0.3924 0.3827 0.000007441
[21,] 0.4297 0.4198 0.000007071
[22,] 0.3921 0.3827 0.000007435
[23,] 0.3657 0.3580 0.000007252
[24,] 0.3654 0.3580 0.000007191
[25,] 0.3654 0.3580 0.000007185
[26,] 0.3402 0.3333 0.000006951
[27,] 0.3527 0.3457 0.000007047
[28,] 0.3279 0.3210 0.000006811
[29,] 0.3402 0.3333 0.000007012
[30,] 0.3417 0.3333 0.000007006
[31,] 0.3400 0.3333 0.000005643
[32,] 0.3017 0.2963 0.000006553
[33,] 0.3219 0.3210 0.000006848
[34,] 0.3149 0.3086 0.000006674
[35,] 0.3651 0.3580 0.000007125
[36,] 0.4040 0.3951 0.000007473
[37,] 0.3118 0.3086 0.000006732
[38,] 0.0480 0.0494 0.000001459
[39,] 0.0130 0.0123 0.000000384
[40,] 0.1076 0.1111 0.000002996
[41,] 0.2192 0.2222 0.000005450
[42,] 0.0361 0.0370 0.000001124
[43,] 0.2162 0.2222 0.000005454
[44,] 0.1540 0.1605 0.000004245
[45,] 0.2311 0.2346 0.000005600
[46,] 0.1905 0.1975 0.000004994
[47,] 0.1557 0.1605 0.000003948
[48,] 0.1912 0.1975 0.000004937
[49,] 0.0584 0.0617 0.000001819
[50,] 0.1782 0.1852 0.000004757
[51,] 0.1916 0.1975 0.000004823
[52,] 0.0126 0.0123 0.000000379
[53,] 0.2365 0.2346 0.000005588
[54,] 0.3652 0.3580 0.000007077
[55,] 0.1946 0.1975 0.000004996
[56,] 0.3270 0.3210 0.000006850
[57,] 0.2378 0.2346 0.000005625
[58,] 0.1572 0.1605 0.000004240
[59,] 0.1969 0.1975 0.000004981
[60,] 0.1682 0.1728 0.000004511
[61,] 0.2484 0.2469 0.000005539
[62,] 0.1970 0.1975 0.000004994
[63,] 0.1101 0.1111 0.000003116
[64,] 0.1975 0.1975 0.000004976
[65,] 0.1240 0.1235 0.000003413
[66,] 0.2221 0.2222 0.000005329
[67,] 0.2347 0.2346 0.000005662
[68,] 0.2092 0.2099 0.000005206
[69,] 0.2205 0.2222 0.000005454
[70,] 0.1845 0.1852 0.000004758
[71,] 0.1722 0.1728 0.000004508
[72,] 0.0486 0.0494 0.000001476
[73,] 0.1474 0.1481 0.000003971
[74,] 0.0844 0.0864 0.000002439
[75,] 0.0843 0.0864 0.000002455
[76,] 0.0721 0.0741 0.000002123
[77,] 0.0478 0.0494 0.000001465
[78,] 0.0240 0.0247 0.000000741
[79,] 0.1798 0.1852 0.000004735
[80,] 0.1187 0.1235 0.000003268
[81,] 0.0480 0.0494 0.000001427
[82,] 0.0454 0.0494 0.000001481
attr(,"gstari")
[1] FALSE
attr(,"call")
localG(x = cpi.df$mo2014m3, listw = nb2listw(nb672km, style = "B"),
return_internals = T, GeoDa = T)
attr(,"class")
[1] "localG"
s2.df <- data.frame( V1=0, cpigetis= 0)
s2.df
for(jj in 1:30) {
nb <- dnearneigh(cpi, d1=0, d2=662+ 10*jj, row.names = IDs)
# Spatial filter for log.GDPPCPCXX or for RGDPPC20XX
cpigetis <- localG(cpi.df$mo2014m3, nb2listw(nb, style="B"), return_internals = T, GeoDa = T)
#the output is created so that it resembles the output of the PPA program
aj <- cbind(cpi.df$IDCity, cpigetis )
aj<- as.data.frame(aj)
aj
n= c(NA, NA)
m= c(jj, NA)
nm<- data.frame( n,m )
names(nm)[1] <- "V1"
names(nm)[2] <- "cpigetis"
aj <- rbind( nm, aj)
s2.df<- rbind(s2.df, aj)
}
s2.df <- s2.df[-1,]
#
head(s2.df)
tail(s2.df)
#
s2.df
locg<- s2.df%>%
select(1,2)
locg
names(locg)[1] <- "point"
names(locg)[2] <- "getis"
locg
in the following code insert the number of points (regions, 514 districs) and the number of different distances used for calculating the local Gi statistic (100 as written int he for loop on line 228 of this code “for(jj in 1:100) {”…)
nump=82
numd=30
rows= (nump+2)*numd
rows
[1] 2520
As you can see this data is very messy, there are two many rows we need two add a new coloumn (r) which refers to the distance used for calulating Gi.
note: just the absolute value of Gi is needed
r<- c(1:rows)
for (n in seq_along(r)) {
for (j in 1:numd){
if ((nump+2)*(j-1)< n & n<= (nump+2)*j ) {
r[n]=662+ 10*j
}
}
}
df <- data.frame(cbind(locg, r))
df<- df %>%
select(point,getis,r) %>%
filter(!is.na(point)) %>%
mutate(getis= sqrt(getis*getis)) #just the absolute value of Gi is needed
df
NA
now having the new column r is possible to spread the dataset so that each observations is the getis statictic each distance is a row and each column represents each point
df <- df %>%
select(point,getis,r) %>%
spread (point,getis)
df
the distance delta is the distance for which the statistic Gi starts to decrease in absolute value
df
dfmat<- as.matrix(df)
dif= matrix(data=1000, nrow = numd-1, ncol = nump+1)
for (i in 1:(numd-10)) {
for (j in 1:(nump+1)) {
if ( dfmat[i+1,j]-dfmat[i,j]< -0.01 & dfmat[i+2,j]-dfmat[i,j]< -0.01 & dfmat[i+3,j]-dfmat[i,j] < -0.01 & dfmat[i+4,j]-dfmat[i,j] < -0.01 ) {
dif[i,j]<- dfmat[i,1]
}
}
}
#deleted this line in the if condition && dfmat[i+2,j]-dfmat[i+1,j] < -0.1
#&& dfmat[i+2,j]-dfmat[i+1,j] <= 0
difframe<-as.data.frame(dif)
difframe
NA
the cells with values equal to 1000 do not have any meaning the other cells represent the distances for which the statistic decreases in absolute value.
for each column there may be many distances for which the statitic decreases in value. Just the first of such values is needed.
the critic value of lambda is at the end of the code (the mode of the distances distribution)
mind<- c(1:nump)
for (i in seq_along(mind)) {
mind[i]= min(dplyr::pull(difframe, i+1)) #Just the first of such values is needed
}
point<- c(1:nump)
hist2<- data.frame(cbind(point, mind))
#hist2
p<- hist2 %>%
filter(mind<1000) %>%
ggplot()+
geom_histogram(mapping = aes(x = mind), binwidth = 10)
p
library(plotly)
ggplotly(p)
hist2 %>%
filter(mind<1000) %>%
summarise(mean= mean(mind), median=median(mind), mode=mfv(mind))
NA
NA
NA
#delta is the mode = 338 km
END
sessionInfo()